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ABSTRACT 


A comparison of standard z transform and algebraic 
substitution synthesis methods for digital filters are 
presented with emphasis on the frequency domain character- 
istics. A generalization of Martin's procedure, which leads 
to recursive filters, is obtained. Also a new method to 
synthesize a digital filter from a magnitude squared versus 


frequency characteristic specification is presented. 
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Lon NP RODUCTIGn 


Digital filters have become increasingly important 
because of technological developments in integrated circuits 
and digital techniques. This progress has increased the 
mumber of applications in which digital filtering techniques 
@an be used. In spite of the fact that it is a Welatively 
new development, a vast body of literature is available for 
digital filters. The relationship >etween digital filters 
and classical continuous filter theory can be depicted as 
mie Fig. lil. 

Start from the differential equation at the center of 
the diagram. The Laplace transform of the differential 
equation may be taken and a transfer function formed. The 
frequency spectrum of the continuous filter may be evaluated 
as shown in the diagram by letting s = jw. The upper part 
of the diagram shows the relationship between the transfer 
function of the continuous filter in the frequency domain, 
and the frequency domain characteristics of the sampled impulse 
mesponsesoL the continuous filter. In deriving these results 
use is made of the Z transform. The lower part of the diagram 
shows numerical integration of the differential equations. 
Adopting a numerical integration algorithm yields a difference 
equation which can generally be expressed as a discrete time 
equation. For a given integration scheme, as illustrated in 
the diagram, it is possible to substitute a function of z 


(depending upon the type of integration) for the variable s 
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of the Laplace transformed continuous equation. When z = 

exp (jwt) is substituted in this equation a spectrum is obtained 
which, in general, 1S not identical to the spectrum of the 
sampled imoulse response of the continuous filter. One of the 
objectives of this thesis is to compare numerical integration 
with the standard Z transform results as depicted in the upper 
Pertion of the diagram. 

The Z transform technique is described in the work of J. 
Ragazzini and L. H. Zadah [1] and has been discussed by many 
other authors including B. C. Kuo [2]. J. Salzer [3] has 
Wereten a Classic paper on the comparison of digital imetegra— 
tors in the frequency domain with the ideal continuous 
integrator. The bilinear substitution, s = = —= , which 
corresponds to trapezoidal integration, was considered very 
early by A. Tustin [4] and has been discussed thoroughly by 
re Steiglitz [5]. 

im this thesis several areas of digital filter theory ane 
considered in detail. 

Chapter two contains a derivation of basic Z transform 
theory, a review of the synthesis of digital filters uSing 
Standard Z transforms, and a discussion of the effect of 
Sampling frequency on the frequency spectrum. The relationship 
between the s and z plane representation of the equations is 
presented. Chapter three contains a general interpretation of 
the algebraic substitution, s = f(z), for several numerical 
integration techniques, a discussion of the bilinear substitu- 
tion (trapezoidal integration), and a development of the effects 


of sampling time on the precision of numerical calculations. 
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Chapter four discusses the synthesis cf digital filters enon 
frequency spectrum characteristics. It contains a derivation, 
discussion, and a practical application of a direct method 

of synthesis originated by Martin [6]. A new generalization 

of this method to obtain recursive filters from a specification 
of magnitude versus frequency characteristics, is presented. 
The advantages of this technique include the fact that a zero 
phase shift filter can be obtained. Derivation of another new 
method Of finding the coefficients of a digital filter, starting 
from a magnitude squared criterion, is presented. The advan- 
tages of this technique include reducing computer storage 
requirements, and reducing the number of multiplications to one 
half of those required by the Martin method. Chapter five 


presents a summary and several questions for future research. 








Li geewlHE 4 TRANSFORM ATE CoherOU ra eee 
TO DISCRETE TIME SIGvAnS ANDSDIGI Eh ares 

Systens may be classified as linear, time invariant, 
causal and so on. They may also be classified according to 
the nature of the signals present. A continuous Signal is 
expressed as a continuous (or piecewise continuous) function 
of the independent variable, t. The signal must be uniquely 
defined at all values of t within a given range, except 
Pessibly ata denumerable set of points, as for example in & 
Square wave. A discrete signal is defined only at a sequence 
of discrete values of the independent variable t. A quantized 
Signal is one which can assume only a denumerable number of 
different values, but quantized sicnals may be either continuous 
=. discrete, as discussed by P. M. DeRusso et al, [7]. 

There are many examples of continuous systems in nature, 
and they are to be found almost everywhere. Examples of 
discrete and quantized systems are perhaps more difficult to 
conceive but this does not mean they are scarce. For instance, 
the output of a continuous system when sampled for digital 
computation, or other purposes, becomes a discrete system. 

Some systems are discrete by nature. The distance to an object 
as measured with a pulsed radar; or a pulse coded communica- 
tions system which is time multiplexed are examples of inherently 
discrete physical systems. A digital computer ts an excellent 
example of a discrete, quantized system. The quantizing levels 


are determined ultimately by the finite number of bits of the 








number representation. The discrete characteristic is deter- 
mined by the fact that the independent variable assumes valves 
at finite increments only. 

In order to deal with discrete time quantized systems, 
existing Z transform theory has been extended [Ref. 2]. This 
theory was developed primarily for sampled data control 
Systems and resembles Laplace transform theory for continuous 
systems in that it iS operational in nature. The Z@ transform 
has proven to be useful for the study of all discrete time 
systems and it is extended approach which is reviewed in this 
Chapter. Several important fundamentals are developed anda 
basic miSsépplication of the theory which occurs commonly in 


the literature, is discussed. 
\ 


. DEP TNT EEON OF THE 42 TRANGYORM OF A DISCRET=S SiGNAL 
Conender a continuous signal x(t), for t>0, when is 
sampled at a frequency, fo or equivalently, every T = = 
seconds. 
Assume that the signal is sampled ideally, that is, the 
Sampled signal assumes the value of the continuous signal with 
infinite accuracy at discrete values of the independent 


miataoles t = nT; n= 0,1,2, «-- 


The process of sampling can be described as 
6 


x* (t) x(t) s(t) (2.1) 


where 


ale 


i 


the sampled signal 
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x(t) = the continuous signal 
s(t) = the sampling function, representing a train of 


impulses, T seconds apart. 


s(t) = = 6 (t-nT) 2.2) 
n=0 

Prom (2-2). and (2.2) 

x*¥(t) = 2 x(n (e—nT) (2.3) 


n=0 


Taking the Laplace transform of the sampled Signal we 


have from Eq. (2.3) 


X*(s) = & {x*(t)]} 
=e), Gs (en (2.4) 
n=0 
Buc £ {§(t-nT) } = een 2.5) 


So that Eq. (2.4) becomes 


X*(s) = © x(nT) e- (2.6) 


Now, define a new variable z, such that 
7 Gabe (2-7) 


so that Eq. (2.6) becomes 


X(z) = X*(s) = YE x(ntT) z” (2.8) 
n= 


sT 0 
Z=e 


The standard procedure to find the Z transform of a 


discrete signal can be summarized as follows: 


ibe 








(a) Sample the continuous signal 

(b) Find the Laplace transform of the sampled function 

(c) Substitute aoe = Z in the Laplace transform of the 
sampled signal to obtain the 2 transform. 

Using common notation, the above can be written as follows: 


ea) OS eee Se [ec ia Gey (2.9) 


sdb 
Z=e 


This process of obtaining the 2 transform of a signal is 
melmstrated in Fig. 2.1. 
aL Obeaining the eazy teas pornueenne  Oosed Born 
The foregoing definition of the 2 transform provides 
a method for obtaining the transform of a Signal, but it has 
two disadvantages: 
a Results are presented in open form because the 
transform of a signal is given as a summation of inverse 
powers of z whose coefficients are the values of the signal at 
the sampling times. For sampled functions derived from contin- 
uous functions with poles in the left half s plane, such 
eer ficients will eventually decay to zero. If the decaying 
process 1s slow compared to the sampling interval, a large 
number of terms in the sampled Z transform expansion are required 
to represent the successive values. 
lo) A The process of obtaining the Z transform may 
be very long, particularly when only the Laplace transform of 
the continuous signal is known. The process requires taking 


the inverse Laplace transform of the continuous signal, sampling 


eZ 








wioyjsuerzay 2 suR uteqqo of sahpssord 


WTOJSURIYL 


soe{TdeyT zetdues 


“tana *Bta 


WIOFSeAL 
soe T dey 


VSTAVAUT 





di 








the resulting time function, and finally expressing the sample 
values as a series of inverse powers of Zz. 

There is an alternate method for obtaining the Z 
transform of a signal that is particularly suitable when the 
Laplace transform of the continuous Signal is known. This 
latter approach gives the result in closed form. From Eq. 


(2.1) 
X*(s) = & {x(t) s(t)} (2.10) 
and so X*¥(s) =x(s) * S(s) 2.11) 


where the asterisk denotes convolution. 


On yee 
X*¥(s) = as | Ke) acs (Sa) ad 2. 1a) 
c-35 


now S(s) = | ME Gotes oWlin) Bae dt 
Oo n=0 
S(s) - > e nts 
=0 
S(s) =1 / (l-exp(-sT)) §2:13) 


provided |exp(-nT) |< 1 
In (2.11), use has been made of the sifting integral 


property of the impulse, namely 


| Ce OME) sete) elie eae (Ocak), (2.14) 
0 
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Therefore, the Laplace transform of the sampled 
function can be written in terms of the complex convolution 


integral 


1 C+} 1 
X* (s) = | x (%) ;  =(saan ar (22.1 5) 
= 


2nj ¢c-jo 





If ry as shown in Fig. 2.2 is taken as the path 
enclosing the left half of the s plane then integral (2.15) 


may be divided into two integrals és follows 


OP di C7 Fo) 


Note that ry 1s composed of the straight line joining 


the points (c,-j~) and (c,+j~%) and the counter-clockwise semi- 


eure le Or infinite radius. 


it wine sCOMGd Le on 


lam s X(s) = 0 @ 17) 


S > © 
is met, then the second integral tends towards zero and Eq. 


(2.16) becomes 


dr (2.18) 


A very important case occurs when the Laplace 


transform X(s) of the continuous signal is a ratio of two 
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The contour integration path 
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two polynomials 


_ Nts) 
X(s) = wT (C28), 
In this case the function, X(s), has a finite number 
of poles and the integral around the closed path can be 
evaluated using the residue theorem. Thus Eq. (2.18) becomes 
SO NO) 1 
X*¥(s) = 2 SS 2 a0) 


n=l D'(A,) 1- oe (S~Ay)T 


where k is the number of Simple poles of X(A) located at dnl 


with 


D'(A_) = ——- G22 1) 


Ped substituting Eq. (2.7) into (2.20), yields Xz) the z 
transform of x(t) 


. N(A)) 


x72 = 


1 . 
; (2522) 


n=1 D' (i) jr a “Rika * 
The above formula (2.22) represents a closed form 
for the Z transform of the signal which contains as many terms 
as the degree of the denominator polynomial of the function 
X(s). Expansion of Eq. (2.22) ina series of inverse powers 


@i Zz yields Eq. (2.8). 


B. SYNTHESIZING A DIGITAL FILTER WITH THE STANDARD Z TRANSFORM 
In the previous section two methods have been developed 
to obtain the Z transform of a sampled signal, starting £rom 


either the continuous time function or its Laplace transform. 


any 








The same procedure can be followed to obtain the Z transform 
corresponding to the sampled output of a continuous filter, 
considering the signal to be the impulse response of the 
meter, h(t). 

in general there are mamy Known procedures fer the design 
of continuous filters to meet various types of specifications. 
These procedures are usually described in the s domain or in 
the frequency domain, s = jw. The Z transform is used to 
translate continuous time realizations into discrete time 
mri zZzatious Or digital filters. 

In control theory, when the differential equations of 
the systems are known, it is often necessary to Simulate 
continuous systems on a digital computer. One approach is to 
consider the impulse response which is the inverse Laplace 
transform of the system's transfer function... The Z transform 
can be applied to this signal after sampling, yielding a discrete 
Simulation of the system. 

Equation (2.8) assures that the Z transform of a digital 
filter or system, when expanded in terms of inverse powers 
of z, represents identically the impulse response of the 
continuous filter or system sampled at intervals T. 

If X(z) is the Z transform of the rnput signal to a system, 
and if H(z) is the 2 transform of the system transfer junction, 


then 


Wezel) HZ) (2223) 
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Here ¥(zZ} represents the Z transtomeort tuewou epi mooie 
and corresponds to the numerical convolution of the input 


Signal an the system transfer function, as is demonstrated 


below — 
7) a x (mT) Bat 
m=0 
7 (2.24) 
fc) ee hc 
k=0 


where h(kT) is the sampled impulse response of the system. 


Subst belting, (2.24) 2mee. (2.2) 


vie) Sor OR sia Tagen = oe (2.25) 
k 
and replacing subscripts, such that ktm = n or k = n-m, then 
the limits of the summations are converted as follows: 
k = 0 becomes n =m and Kk = © becomes n = ©... 


Thus (2.25) becomes 


Cc Cc 


Y(z) =X © xl(mt) h((n-m)T) z™ (eae) 


m=O n=m 


Changing the order of the summations and neglecting or 


adding zero terms of the type h(-kT) (k=1,2,.....) yields 


co n 
Y(z) =2 5 X(mT) h((n-m)T) z”™ (Doan 
n=O m=0 
= y(nT) z 
n=0 


1g, 








Hence the coefficient y(nT) of Zz in Y(Z) is given by 
n 


y (nd) = 2 <x Gia (4) (2283 
m=0 


Equation (2.28) is the numerical convolution of the input 
samples and the digital filter sampled impulse response. 

The output signal transform, Y(z), when expanded rs pipes) 
series of inverse powers of z represents the output signal at 
different times. If X(z) and H(z) are obtained from sampling 
a continuous input Signal and the impulse response of a con- 
tinuous system or filter respectively, then the coefficients 
in the expansion of Y(z) will correspond to the sampled output 
of the continuous system or filter. Hence when synthesizing 
amaigital filter using the Z transform of a continuous filter 
impulse response sampled at intervals T, the resulting output 
Signal is an exact description in the time domain in the sense 
that the output values of the digital realization will be equal 
to the values of the continuous filter output signal sampled 
at intervals T. 

Regarding the frequency spectrum of the filters, consider 
pme@eontinuous system with transfer function H(s). If an input 
Signal E, (s) 1s applied, the output signal E. (s) has a Laplace 


transform 
E, (s) = H(s) E. (s) (2229) 


The inverse Laplace transform of E,{s) is the output 
Signal e, (t) = £-" {ze (s)}, which now is sampled at the rate of 


1 / T samples per second, yielding the sampled output signal 
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e,* (t) as shown in Fig. 2.3a. The Fourier transform of this 


Signal is the desired output spectrum 
* (5 = *« 
E * (ju) = 9 fe *(t)} (2..30) 


This spectrum could also have been obtained by @Genvolvang 
the Fourier transform of the sampling function ys 6f7 (a = 
n=-o 


nw.) ] Suechethat 


co 


E,* (jw) = EB, (jw) = > 6 [5 (w-nw .) ] (2%,.3 1) 
= J EB, (j fw-nw ]) 
n=-© 


This signal processing may be compared with the develop- 
Mmemt Of the digital filter shown in Fig. 2.3b. We assume Ene 
‘same sampling frequency f = 1 /T and obtain the Z transform 
faez) Of the filter function H(s). Then the Z transform Gf Ehe 


meter output E, (2) is given by 


E,(z) = H(z) E, (z) | (2.32) 


where E, (z) is the Z transform of the sampled filter input 
e,*(t). Substituting 2 = exp(jwT) into E, (2) yields the desired 
output spectrum E,* (34) - The following question may now be 
asked: Is this output spectrum of the digital filter identical 
wen the output spectrum obtained in Eq. (2.31)? This questron 
is discussed in the next section. 

lie The Effect of Sampling on the Prequency screceamm 

For a signal of finite bandwidth w radians called a 


band limited signal shown in Fig. 2.4a Shannon's Sampling 


72) 
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Theorem [8], requires a sampling rate equal to at least twice 
the maximum frequency w of the signal in order to recover the 
Signal without distortition after passing it through an ideal 
low pass filter. The spectrum of the signal sampled ata 
higher rate than the Nyquist rate is shown in Fig. 2.4b. If 
the sampling frequency is lower than the limit set by Nyquist, 
then some overlapping of spectrum will occur. This is known 
as the aliasing effect and is illustrated in. Hg. 2.4c. It has 
been shown by L. J. Fogel [9], that -it is possible to recover 
a Signal sampled at a lower rate than Nyquist (or one for which 
aliasing has taken place) if for every Sample additional infor- 
mation about the signal is Known. For example, this is possible 
if its derivatives with respect to time are known. However, 
this idea cannot be applied here because we restrict ourselves 
to single input filters where only the values of the signal are 
known. 

To answer the previous question regarding comparison 
of the two spectra of Fig. 2.3, consider the continuous case 


of Fig. 2.3a, where for s = jw 
BE, (34) = H(jw) E, (jw) (2.33) 


After sampling this spectrum, E,* (4) is the convolu- 
Elon Of EB, (34) and the spectrum of the sampling signal, which 
is a train of impulses separated by w. on the w axis as given 
by Eq. (2.31). The effect is to make the sampled spectrum 
repetitive with period Wo From Shannon's Sampling Theorem, 
the sampled output can be reconstructed without aliasing distor- 


tion provided E (5) is band limited and the sampling frequency 
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(a) 


(b) 


(e) 


ElGiar 2a4% 


JE (ju) | 


Spectrum of a band limited signal 


JE (jw) | 


Spectrum of a sampled band limited 
Signal, fia > 2w 





aliasing 
effect 


Spectrum of a sampled band limited 
Signal, 2nf < 2w 


Effect of sampling frequency on the 
frequency spectrum of the sampled 
Signal. 
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is higher than two times the maximum frequency of Eo (jw). For 
E, (5) to be band limited either E, (jw) or H(jw) have to be 
band limited. 

Consider now the discrete case shown in Fig. 2.3b, 
where e, (t) 1s sampled and the Z transform is obtained. From 


mee 9 (2.8) 


E, (2) == E.*(s) (2.34) 


Substituting Eq. (2.34) inte @equation (2.32)) yietds 
Bo*(el@*) = weet") BL *(e7%*) (2.35) 


JwT 


As may be seen, the outptt spectrum E,* (e ) 1s the 


product of the two spectra H* (et) ana E.*(e7"). Clearly 


both these spectra must be undistorted (no aliasing) if 
Eo* (eI) is to be undistorted. Hence for this case the 
requirement is different from the case of Fig. 2.3a. In view of 
the above remarks we conclude that it 1s necessary to have 

E, (jw) and H(jw) band limited and furthermore, one must sample 
at a frequency above the Nyquist rate for both spectra of 

Bags. 2.3a and 2.3b to be identical. 

As an example consider the design of a digital low 
pass filter with a cut off frequency of 2 KHz which is required 
Memiilter a signal of 10 kHz bandwidth. Such a filter must 
bee adesigned with a sampling rate greater than 20 KHz. Nete 
that if the signal had been obtained from a continuous filter, 
it could be recovered by sampling at a rate of only 4 kHz. 


It is clear from the above remarks and the example given, 


that there is a maximum frequency that can be handled by a 


Zo 








digital filter in real fEime operation. | 0in, face thi se mins 
depend on the computing speed of the processor and the general 
eomplexity of the filter. It is destnable to ccoo sig mea 
miter designs as simple as possible, although) comamdaratiom. 
of stability may modify these considerations. 
Jl The Interelationship Between the s and z Planes 

The interelationship between the s plane and the z 

plane can be derived from the definition of the variable z 


namely 
z=e C2ie Sie) 
A point in the s plane, namely (o+jw) will he 


mapped into the point in the z plane 


— o (otjw)T x Qo _JuT (2.37) 


Therefore, the points in the s plane with constant 


camping, 1.@., 6 = constant, map onto a circle of constant 
radius, ae in the z plane. Points in the s plane with 
Constant frequency, i.e. w = constant, map onto radial lines, 
arg 2 = wT. 


Furthermore, for points in the left half s plane 


o < 0, and so 


2B) = el <1 for T> 0 (2.38) 


Equation (2.38) indicates that the entire L.H.P. 


Maps into the inside of the unit circle in the z plane. 
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Although useful, the foregoing discussion may be 
misleading unless it is applied proverly. The misuse of 
He, (2.36) occurs often in the literacune.) Take forme amos 
mae paper by C. J. Greaves and J. 7A. Gadze. [10]. Taeietidie= 
take consists of replacing the variable s in the Laplace 
ieamsform of the continuous filter by 2s corresponding Z, with 


s = (1/T) ln z as 


GZ ee eS) (2...39a) 





So- G ty sin -Z 


This eguation is in error and should be replaced by 


HAZ) = HX (Ss) (2 53 2D) 





s = (1/T) ln z 


The reason is clear from the closed form expression 
for the Z transform, Eq. (2.18), where before the substitution 
of variables is made, a complex integral has to be evaluated. 

However, there is a direct connection between the 
Laplace transform of the continuous function and the Z trans- 
Zea or the sampled function. Namely “the poles of the Laplace 
transform of the continuous function are mapped into the z plane 
of the sampled function using the mapping function of (2.36). 

This can be seen from equation (2.22) where every 


term in the summation will have a pole in the z plane at 


Z=e (2.39c) 


which indicates the mapping of the singularity es accomding to 


the mapping relation of (2.36). 
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In summary Eq. (2.36) may only be used to transform 
poles of the continuous function in the s plane to pole 
locations of the sampled function in the z plane. It does not 
mepresent a chenge in variable, except swnen app litedma em alin, 
to a sampled signal. As noted this subtlety has proven to he 


a source of many mistatements in the literature. 
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Lil. THE ALGEBRAIC, SUBSTEULION MEtHoObD 


An equivalent digital representation of a continuous 
system Or filter is often required in practice. Tharescould 
arise, for example, in a digital simulation of a continuous 
control problem, or in an attempt to extend continuous filter 
theory to provide digital filter realizations. This latter 
application is very appealing because continuous filter theory 
is very well developed and aeeedhanes exist to match given 
specifications to filter relizations. Designs have been 
tabulated and are available in many handbooks. If a direct 
method to obtain the equivalent digital realization of a contin- 

/ 
uous filter exists, it would make all previous knowledge of 
‘continuous filter theory available for the discrete case.. 

It 1S convenient when implementing a digital filter, to 
obtain its representation as a Z transform, and particularly 
as a ratio of two polynomials of inverse powers of z. Such a 
representation indicates the linear combinations which have to 
be performed upon past values of input and output signals 
samples, to obtain the present value of the output. 

ime most common description Of Continuous falters Tipehews 
domain iS aS a ratio of two polynomials in s and so it is con- 
venient to have a rational function of 2 to directly replace 


for s in the continuous type description. This is if 


ee) (22a) 


then H(s) = H(f(z)) = G(z) . (3.2) 
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If the fumetiongi(Z).1s.- a lratlo.Orestwe spol nemaucscme:. 
inverse powers of z and if H(sS) is a ratio of two polynomials 
in s, then G(z) in Eq. (3.2) iS assured to be the ratio of two 
Polynomials an inverse powers of 9Z- 

This procedure, described in Eq. (3.2) is called the 
algebraic substitution technigue, and converts rational 
polynomials in s into the desired rational polynomials in z. 

It will be shown that the selection of a specific f(z) 
to substitute for s, corresponds to the adoption of a numerical 
integration algorithm to integrate the phase state variables 


of the filter or system in discrete time. 


A. MEANING OF THE ALGEBRAIC SUBSTITUTION s = £(z) 

In order to find functions that can replace the variable 
fom the description of the continuous filter, H(sjigethe 
meaning of the substitution must be understood. In the first 
place, the variable s of the Laplace transform theory represents 
ame operator in the time domain, i.e., it is a differentiator. 
Similarly, the inverse of s, a. represents an integrator. 


Consider a filter whose transfer function is 


cpecineel, 8 ritcicient 5” 
H(s) = ———————_______ , where m > n (3x) 
beet best eee bes 
O ik m 


If this transfer function is the ratio of the output signal 
transform to the input signal transform with zero initial 


conditions, then Eq. (3.3) can be expanded as follows 


Y Vis) acs eens) (3.4) 
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The following equations can be formed arbitrarily 





W(s) _ i 

8 (3.5) 
X (s) b + by s + b. s+ 4 b S 

Woks eh ae Ti : 
ATES Woe ao + a, s + ... + a, s (3:.6) 


Rearranging Eq. (3.5) and (3.6) yields 


m-L . 
b s" w(s) = X(s) - 5 b. s? W(s) iS. 7) 
m ‘ 
jee 
nN ; 
Y(s) = = a.s? W(s) (3.8) 
j=0 


EGQuations (3.7) and (3.8) can be aliustrated as angie: 
3.1, indicating a possible implemeritation, as is usually done 
‘in an analog computer simulation for low order systems. 

Tt has been shown that a transfer function of the type 
described in Eq. (3.3), can always be expanded in the form shown 
maerig. 3.1. From this illustratien it can be seen, that the 
replacement of the integrations (denoted by 1/s in the state 
ee expansion of a Continuous E£ili@ermim by a function of Z, 
1/f£(z), corresponds to a direct algebraic substitution of the 
variable s by the function s = f(z) in the closed form expression 
Srebq. (3.3). 

Furthermore, the linear operation I/f(z), should perform 
an approximate operation in the discrete time domain similar 
to the operation, l/s, of the transform domain, if the responses 
obtained with the digital filter are to resemble the responses 


obtained from the continuous system. 


Sid. 
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In view of the. above remarks, it can be concluded that 


the algebraic substitution technigue represents the adoption 


on a numerical 


Om 1ntegqracion 


B. OBTAINING 


To derive 


integration algorithm to perform the operaticn 


for the state variables of the filter. 


RECURSIVE INTEGRATION ALGORITHMS 


recursive integration algorithms, it must be 


decided whether or not the present value of the input is 


available to compute the present value of the output. This 


question can be settled if some assumptions are made. For 


iastance, 1f the time required to perform the computations is 


considered negligible, that is if the time elapsed between the 


present input and the corresponding output of the filter is of 


no importance, 


then the algorithms can be derived using inter- 


‘polation formulas. On the other hand, when this factor is 


considered important and the output. at the present time is 


computed without the present input value, then some extrapola- 


tion formula must be used. This latter condition involves some 


form of a predicting method, and intuitively seems to be less 


accurate. Among several different criteria to find recursive 


integration algorithms, the following methods are described 


because they lead to well known integration formulas which are 


Mseed in practice. 


a. The Adams-Moulton Method. 


This method is derived from the expression 


sr. 
ii | ze) cle (3.9) 
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where f(t) LSsja polynomval yor (n=) 75 degree approximating 
the input signal, based on the last n available input values. 
The recursive integration algorithms fom thes os 


two degrees of approximating polynomials are 


“Interpolating Form Extrapolating Form 
fee tUsing Last input) (Not using last input) 
2 a = Ty = 
meee k= 'D x tXy-)) hee Se 7, See! 


“a Tl oecy _ | = ms = : 
Be ka) 1a ga ee? Yeo Oe OX 


From these algorithms it can be seen that the first 
degree approximation uSing the last input value is the well 
mewn Trapezoidal rule of integration. 

b. The Milne's Method. 
This method correspond to the numerical integration 


Gri cerion 


+ | f(t) dt (3 GF 
k=-n 


where f(t) is the polynomial of gees 


degree approximating 
the input function using the last n available values of the 
inputs. 

Clearly, when the approximating degree n is unity 
Adams-Moulton and Milne's methods are the same, giving identical 


recurSive integration algorithms. For a polynomial approxima- 


tion of second degree the recursive Milne's algorithms are: 


>» 


ore ee 
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Interpolating Form EXtrapolating Forum 


n (USing last input) (Not using last input) 
2 Yr1=yY fe (X, +4X = ) Y,=yY +2! (ox -X +2X ) 
k “k-2 3° k sk eS kK ke se Nea dl c = ke3e 


It can be seen that Milne's first degree formula also leads to 
the trapezoidal rule and that the second degree formula leads 

to the Simpson's 1/3 rule of integration, when the last input 

is used. 

These recurSive integration algorithms as well as others 
-have been analyzed in the literature under the general category 
@£ digital integrators, as for example in one of the first 
papers on the subject which was published by John M. Salzer 
mmr 953 [3], or in a paper by A. P. Sage and R. W. Burt Jia]: 
| To obtain the frequency spectra of the algorithms, one 
must first obtain the equivalent Z transform representation. 


For example consider the trapezoidal rule of integration 
(3.11) 


The equivalent digital filter can be derived by taking 


the Z transform 





aim) = a = Sle) & t (X(z) + ab eee) (3.12) 
and the filter transfer function is 
ViGZ T 1 + zt 
Haz) = = — i (3.23) 
X (Zz) 2 1-2 
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This function corresponds to the algebraic substitution 


= 43.14) 


The frequency spectrum of the filter described in ((3.13) 


is obtained by substituting z = exp (jwT) 


; 50 
(Clue), oo (3 Es) 

3 ON eee 
a1 (Qoee -3 5 cotan = (3.16) 

and hence 
bH(e2 OF) | = cotan oS Re 17) 
jwT = TT G 

Arg H(e )=-+-r 3.18) 


These equations can be compared with the magnitude and 


phase of an ideal integrator 


His) = = (3.19) 
H(jw) = = (3.20) 
jH(jw)| = ~ (3.21) 
Arg H(jw) = =4 (22 2)) 


Equation (3.17) can be expanded in a series of the form 


SI nb ee - (3.23) 


ey ag Z wT wp J, wr 
2 WT 6 360 on W 1 
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and therefore the magnitude spectrum of the trapezoidal rule 
digital integrator approximates the magnitude of an ideal 


mabegrateor for yalueseor Ww Suchmenas 





<< = (3.24) 


W << vie Ce... 25) 
LT 
e1g 
eee ay conEn 
Tt. Ss 2s 


In Salzer's paper a study of the frequency spectra of 
several different digital integrators 1s considered with a 
comparison to the spectrum of an ideal integrator. 

Figure 3.2 shows es of the trapezoidal rule and 
Simpson's 1/3 rd rule plotted against Q = W/W The straight 
line or real axis represents the ideal integrator and as can 
be seen the trapezoidal rule falls below this line, while 
Simpson's rule lies above this line. 

In general it can be said that the more complex the 
integration algorithm, the wider the band of coincidence with 
the ideal integrator and also the larger the number of multi- 
plications to be performed or the increase in computation 
time. As a result a possible figure of merit might be the ratio 
of 3 db bandwidth to the number of multiplications. As a matter 
of fact, simpler integration schemes perform as well or better 


than the complex ones in terms of the above figure of merit. 
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Cc. THE BELENBAR cube Teuton 
This particular case of the algebraic substitution 
technique corresponding to the adoption of the trapezoidal rule 
of integration has a different interpretation p~sds et oeone mero 
Ree Steiglitz [5], and also discussed by Gold and Rader [15]: 
In order to avoid confusion between variables, the 


following notation will be used: 

S=o + jw 29) 
= o 69 ee) | (aeZ.8 ) 
When the bilinear transformation 1s made, such that 


Jae) ar 7 (3.29) 





then the entire jw axis is mapped onto the unit circle itn the 
Z plane. This can be shown from the inverse relationship of 


bee (3.29), namely 


ee Pas 
eee ye eee a (3.30) 
1-(%)s 1- (S)5e 
2 Ded 
To find the relationship between the variable w in the 


s plane and the variable no in the z plane when s = jw, the 


following equation may be used 


T 
ae , -1 of 
a i+4* 5 Jeb . 22 tan ; (3.31) 
Sars) pe 


Equation (3.31) may be reduced to the form 


.- 2 wf Gmse 
fu GUS Sieun 5 Sere) 
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Equation (3.32) indicates Ene transtormation relacraqeene 
frequency axis of the s and z planes. Equation (3.29) means 
that when the bilinear transformation is applied to a rational 
expression in powers of s, H(s), a rational expression G(Z) in 


powers of z is obtained. 


Y - i 
i + 1) (3.33) 


Giz) = Hf 
Furthermore, when the frequency spectrum of the s domain 
representation is obtained by evaluating the expression H(s) 
when s = jw, and H(jw) is obtained, and similarly if the fre- 


quency spectrum of the expression iin powers of z is evaluated 


along the unit circle, the function obtained is 
Beene) = H(Z arc tan >) (3.34) 


The resulting function has the same spectrum as H(jw) but 
eemmressed into the interval -1/T < * < q/Tf. 

The nonlinearity of the frequency scale can be overcome 
by prewarping [Ref. 15] or frequency scaling the chosen contin- 
uous realization, such that after the transformation the 
selected frequency is at the desired value. However, with this 
prewarping technique only one frequency may be scaled properly 
and thus certain characteristics of continuous filters as roll- 
off per octave are not preserved. This method of synthesis 1s 
particularly well suited to synthesize low pass filters with a 
cut off frequency much smaller than the sampling frequency, 


because then the transformation function for the frequencies, 


Ww == tan = may be considered linear, preserving then the 


characteristics of the continuous filter realization from which the 
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characteristics were obtained. This situation can be forced 
by increasing the sampling frequency or equivalently by 
decreasing the sampling interval T, and this solution could be 
the cure for all problems, but this is not so. As was 
discussed in chapter two, reducing the sampling interval, is 
equivalent to reducing the time available for computations, 
and for a fixed computing speed, this results in a reduction 
in the maximum complexity of the filters that can be implemented. 
Therefore, the reduction in the sampling interval is the least 
desirable tool to be used. 
Besides this consideration of computing time, there is 
also another consideration about the number of bits required 
in the digital representation of the coefficients of the filter. 
First it must be noted that the coefficients in the digital 
filters have a finite accuracy, and hence the location of the 
poles in the z plane is made on a discrete basis. To clarify 
this consider a simple example where the location of a pole is 


eamectly related to the coefficient accuracy. The filter 


1 ie 
H(z) = —— = "lias | (Seo) 
A C St Ci Z 
has a pole in the z plane located at the point z = Ci: Because 
of the finite accuracy in the representation of C, this pole 


1 


can only be located at a finite number of positions. That is 
if an bit binary representation of the coefficients is used, 
then the maximum accuracy in the pole position due to quantizing 


effect is 2". 
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Consider then the effect of a seductiensin enc (eamoneena 
interval on the mapping of the s plane into the z plane under 
the bilinear transformation. 

To visualize this, a grid in the s plane, shown in Fig. 
3.3, is mapped into the z plane using the bilinear transforma- 
mae fOr Values of T = 0.1, 0.05, 0.025, 0.1255 The sesuice 
eeee Shown in Figs. 3.4, 3.5, 3.6 and 3.7. 

As can be seen, the area in the z plane into which, a 
given area of the s plane is mapped decreases with T, the 
sampling interval. Conversely a given area in the z plane 
represents a larger area in the s plane as T decreases. 

Consider a point Zy in the z plane distant one quantizing 
level from the point z = 1 (towards the origin on the real 
fers) and distant one quantizing level along a vertical direceven 


(above the real axis) as shown in Fig. 3.8. 
2, = (l - E,) oe En (3.36) 


The corresponding point in the s plane is 


a. a cil ae z ent) bn (see) 
1° Tz tl” T 3-2 45 E, : 


Zoe ee Gt VAG) tea | tee) 
- _ _h 5 (3.38) 
T ak 
vs 
2°) tee) (Amey cae). & 
ae oS RD a seen eae ea Ee 
T (2-F )° + &L (2-F ) + GF 
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Mavping of the grid of the s plane 
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uSing the bilinear trnasformation 
for T = 0.1 sec. 
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Mapping of the grid of the s plane 
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using the bilinear transformation 
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Ze (ig ely) 5 lag 2 


ee: $$$ + jo ——___>—__ (3.40) 
T (2-€)* + €% (2-6,)* + &* 

~ on , Bn : : 

See apy dle ca where En se ee OS (3.41) 


The above derivation can be considered as the mapping cf 
mmesaquare of the Z plane, each side measuring one quantizing 
level, Ent into an area in the s plane of side E,/t: Because 
Sesquantizing effects a singularity at Z., can only be located 
at a finite number of positions, BAen of which corresponds t:o 
a position Si. 

Equation (3.41) indicates that the grainniness of allowed 
locations in the s plane is dependent on the quantizing level 
and the sampling interval. When uSing a continuous filter 
‘design to synthesize a digital filter, the pole locations so 
obtained may not be conincident withthese allowed positions, 
and so it may be necessary to shift them to those allowed 
locations which are closest to the design Locations. Due =e. 
other considerations, like pole location sensitivity in the s 
plane, the maximum shifting allowable in the location of a vole 
in the s plane may be limited to a maximum distance A. If it 
is desired to make the maximum allowable shifting comparable 
to the grainniness in the s plane the following equation may 


be established 


— _ | Yes a 
A 3 (3..42a) 
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In the Z-plane this grainniness would require n bits where 


=n (3.42b) 
From (34.2) 
queens Jere) (saa) (3.43) 


where n represents the number of bits used in the binary 
representation of coefficients and A iS an indication of the 
grainniness of the grid in the s plane. 

The lower bound of the required number of bits can be 
determined from Eq. (3.43). Actually, the required number 
may be more, especially if the poles are located very far from 
the origin in the s plane, because in this case the distortion 
in mapping of areas is bigger and can no longer be neglected 
iieecoemputing A through this derivation. 

From the preceding discussion, it may be concluded that 
the finite representation of coefficients introduces an error 
in the location of the singularities of the function under- 
going the transformation, and that this error increases, that 
is represents a larger change in the s plane when the sampling 
interval is decreased. In other words, it is not possible to 
obtain a better approximation by decreasing only the sampling 
interval. Instead the number of bits which represent the 
accuracy of the coefficients must simultaneously be increased 
to obtain a better approximation to the analog specification. 

Because for real time operation all computations must be 


performed within a sampling interval, the above considerations 
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set a more rigorous limit than previously considered. This is 
so because an increase in bits to represent the coefficients 


tends to increase the computation time. 


D. CONCLUSIONS 

The algebraic substitution method provides an easy method 
to obtain a digital filter which approximates the character- 
istics of a continuous filter, providing then the desired link 
between the continuous filter theory and the discrete filter 
theory, with some restrictions. 

Of all the substitutions, the bilinear transformation is 
the easiest to apply and does not increase the complexity of 
mime Lilter, 1.e€., a ace order continuous Ellter iS transtonrmee 
MEO a ne epaerm ciserere filter. 

The simulation of continuous systems obtained by the 
substitution technique in general can be improved by decreasing 
the sampling interval. 

However, the spectrum obtained with a digital filter using 
this technique is only an approximation to the spectrum of the 
Bees nuous fl beer Evomewniwen 1C isederivea, and 16 themcon. 
tinuous filter is itself an approximation to a given specifica- 
feonm there 1S no assurance that the digital filter will also 


be within the specifications. 


S)8 








IV. sgeYNTHESIS OF DICTITARRE Lie eS 
PROM FREQUENCY SPECTRUM Cla AGT Eo vite 

The standard Z transform technique and the algebraic 
substitution method considered in chapters two and three 
mopnesent an attempt to extend the continuous filter theory to 
discrete or digital filter theory. Gee 1s no method 
available which will guarantee that a digital filter will have 
an exact replica of the frequency spectrum of the continuous 
filter from which it is derived, except when the input signal 
and the filter transfer function are band limited and sampled 
according to Shannon's theorem using the Z transform technique. 
When using the above methods to approximate the frequency 
Peamacteristics of the continuous filters (which themselves 
are an approximation to certain given specifications) there is 
no assurance that the resulting digital filter will meet design 
specifications also. 

A more direct approach to the study of the frequency 
spectra of digital filters seems to be necessary in order to 
provide better answers than the techniques derived from 
continuous filter theory. 

A part of the material presented here is an extension 
of an idea mentioned by M. Martin [6], in which a nonrecursive 
filter is constructed. The approximation of Martin to a given 
specification using a Fourier series truncated at the sien 
harmonic, requires 2 n words of storage and 2 n multiplications. 


The contribution described in this chapter leads to a recursive 


=z 





filter, constructed from a ratio Of wtwo finite Poumiewescmaes 
approximating a given specification. This new procedure 
requires the same storage and number of multiplications, but 
allows one to obtain phase versus frequency characteristics 
other than the linear phase versus frequency with slope obtained 
when uSing Martin's method. 

The last part of the chapter presents che development ores 
generalized technique to obtain recursive digital filters 
derived from a magnitude squared criterion. The magnitude 
squared function 1S expanded as a ratio of two finite Purier 
series. For the approximating series to involve up to the aoe 
harmonic, this procedure requires only n words of storage, 
which represents a substantial improvement over the magnitude 
‘methods discussed in the early part of the chapter. 

A. SrNtE Eos OF NONRECURSIVE FILTERS FROM Avy GIVEN MACNITEZ: 

SPECTRUM 

In certain applications only the magnitude spectrum of the 
output signal is important and a time delay of some sampling 
intervals can be tolerated. The procedure that follows allows 
the given magnitude function to be approximated with a finite 
summation of terms of a Fourier series. This procedure is 
mentioned by M. Martin [6] with reference to data transmission 
filters. 


Consider the elementary filter whose transfer function is 


oe the form 


H(z) = 5 [14277] (4.1) 


A flow diagram of the filter is shown in Fig. 4-la.. 
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(b) The elementary filter of a order 


Fig. 4.1. Elementary filters used for the 
Fourier series approximating technique 
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The elementary filter described has a spectrum of the 


form 
ee = stl + a Jeet 
ere. = exp(jwT) has been 
(4.2) can be reduced to the 
H(eJ4T, = cos wT a ee 
amc hence baie =—wlicos 
Arg H{eJ®*, = - wT 


substituted in Eq. 


£OiM 


wT 


The elementary filter described in (4.1) 


an is 


(4.2) 


EGuatilen 


(4.3) 


(4.4) 


(4.5) 


LS Cenctaered 


to be of crder one and thus a family of elementary filters of 


megner order can be generated in the form of 


{1 + a is 


hol Fe 


H (2) = 


Such filters have frequency spectra of the form 


H (22 *) = cos nwT e 
JH (eI@F) | = |cos nuwT 
Arg aes = - nuwT 


=cut 


(4.6) 


(4.7) 


(4.8) 


(4.9) 


The filter of order n is illustrated in Fig. 4.2b.. 


Using these elementary filters it is possible to approxi- 


mate any given magnitude versus frequency characteristic 


feel" 


Le aS will now be shown, 


First it must be noted that 1f the output of any of these 


filters described in Eq. 


215) 


(4.6) is delayed m sampling intervals, 








the frequency response becomes sor cere rn 


H, (2) = * [l + yey ee as (4 eG 


H(eI°*) EES iene OO es 


—j (n+m) wT 


tl 


cos nwT e 2 


Comparing Eq. (4.11) with Eq) (4.7), 26 is Clear yenoce ene 
Magnitude of the frequency spectrum does not change by delaying 
the output. Only its phase changes in the form of an increase 
in the slope of the linear phase versus frequency 
eharacteristic. 

Therefore, by delaying the out:put signal of these filters 
a convenient number of sampling intervals, it is possible to 
‘make their output signals have the same phase characteristics, 
and hence if several of these specially delayed filters are 
added in parallel, the resulting magnitude spectrum is the sum 
of the individual magnitude spectra because all the complex 
quantities to be added have the same phase. 

An example of the above paragraph is illustrated in Fig. 
4.2a where the elementary filters of zero and first order are 
delayed such that their phase versus freauency characteristic 
become equal to the second order elementary filter. In Fig. 
4.2b the canonical form of the filter is shown, indicating 
symmetry in the coefficients. 

Generalizing, a filter synthesized as a summation of the 


delayed elementary filters in the form 


26 
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(a) Equation form 





(b) Canonical form 
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Synthesis of nonrecursive filter 


approximating a given spectrum function 
involving up to the second harmonic of 


the Fourier Series. 


Si 








(4.12) 


has a frequency spectrum which can be obtained by substituting 


= exp(jwT) into Eq. (4.12), yielding 


: N ' 
Gel) 2 ye cae an Go! (4.13) 

=0 
oT N 

and hence 1G (e4 = dy C cos nuT (4.14) 

n=0 
Arg ele! aa) =e Niwa 4-15) 
For real coefficients, Cur Wall lion tO) eee) oc 


Equation (4.14) indicates the possibility of approximating 
‘any given magnitude spectrum with a Fourier series of finite 
numbers of terms. It is known that the magnitude spectra are 
even functions, and repetitive with period Oo = 21/7 Le. bors 
characteristic assures that it will always be possible to 
expand a magnitude spectrum function as a Fourier series with 
cosine terms only, assuring also the generality of Eq. (4.14). 
Equation (4.15) indicates that the resulting function 
will be delayed N sampling intervals, N being equal to the 
number of harmonics considered in the Fourier series expansion. 
If ‘ Garces nut 20  foreaitas (4.16) 
n=0 ” _ 
then the phase given by Eq. (4.15) will not display any dis- 
continuity jumps and will remain linear. This can always 


be accomplished by proper adjustement of the coefficient Co: 


3ys, 





It is well known that the Gibb's phenomenon occurs near 
a discontinuity of a function. In other words no matteric. 
many harmonics are considered in the Fourier series approxima- 
tion, the Summation of terms does not converge to the value 
of the function near a discontinuity, but does converge to an 
overshoot and a undershoot of approximately 10% of the 
discontinuous jump. In order to make certain that (4.16) is 
Satisfied, the value of Co must be adjusted to include the 
undershoot of the Gibb's phenomenon. In the case of an ideal 
low pass filter the Gibb's phenomenon is about 10% of the 
meesccontinuity. If Co is adjusted so that the Fourier series 
is always positive, this implies that the attenuation in the 
stop band iS approximately 20 db. 

Some other methods may be used to obtain a greater 
attenuation. For instance the Lanczos and Fejer methods to 
attenuate or eliminate the overshoot and ripple of the Fourier 
series can be used as discussed by R. W. Hamming [12]. 

A different approach to get around this Teer ences is 
to work numerically, with the given magnitude spectrum as a 
set of data points, through which a continuous function is to 
be fitted. The lack of discontinuities assures that Gibb's 
phenomenon will not occur. 

A simple example of nonrecursive digital filter synthesis 
is now presented. 

Synthesize a nonrecursive digital low pass filter approx- 
imating the magnitude function shown in Fig. 4.3, where wo = /3T 


and the order of the filter is chosen such that the approximation 


involves the fourth harmonic in the Fourier series expansion. 


a1) 





—— 





“ZOUTTF Te3tHtp e& yyTM 
poyo jeu eq OF OTASTAeROeAeYO ADuSnbezzZ snszsA spnyztubep 


us eS é. 
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Representing the magnitude funcetoneas a Zourgeouesorle] 








yields 
JwTt Ge “ 
|G (e ) | = » Ce Cos .uut 
n 
. n=1 
WW) W) 
where C = a | = COs NWT cu.— 4 i sin nwT | = 
n Ww Cee 
S O S 0 
4 4) sin nTw 
- nTw Sng wie a W nTw 
For this example W/W 5 = 1/6 and wT = 17 os 
Therefore the coefficients Ch become 
Cc = Ay ST Sy ee 
n 3 M173 
Cy = 0.66666 Cy = 10) GSS) ies Cc. = 9.27567 C4 = 0.00 Cy = —~O. 237 ce 


Consider first the case where C is not adjusted so that 


(eo) is not satisfied. 


Jot) | = 0.33333 + 0.55133 cos wT + 0.27567 cos 2uT 


IG(e 


+ 0.0 cos 3wT - 0.13783 cos 4wT 


and from Eqs. (4.14) and (4.12) 


Bij 0.33333 2°) ee 
te atyen? ~ 9223783 og 274 


2 
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and rearranging terms the filter becomes 


2 3 4 


G(z) = - 0.068915 + 0.13783 26 = @227567 2. etOncoee omen 


5 6 8 


+ 0:27567 2 ~~. 0.13783 Ze 20 oGeolo eee 


In Figs. 4.4 and 4.5 the magnitude and phase versus 
frequency characteristics are shown for the filter of the 
example, where the coefficients are obtained directly from the 
Fourier series expansion. The discontinuities in the phase 
function can be seen to correspond to every cross over point 
shown in the magnitude function by an increase in attenuation. 
The lobes in the attenuation band of the magnitude function 
are the ripples of the Fourier series but distorted by the 


Meqarithmic scale. If Co iS ed jusved so@enat Co = 0.46030 then 


2 3 4 


Giz) = = 02068915 + 0.13783 2 ~ £90227567 2 ~ + O-4euccum 


5 6 8 


2 10.27567 zo £eOnlageo) emo OG eo lo 


The magnitude and phase response of this filter are shown 
Maeraigs. 4.6 and 4.7. 

In general the foregoing procedure requires 2 n words 
of storage for n harmonics considered in the Fourier series 
approximation. Furthermore the output signal has a delay of 
n sampling intervals due to the remarkable characteristic of 
Enis type of filters of having a linear phase with a slope 
dependent only in the number of delays of the filter, Failure 
to observe the condition stated in Eq. (4.16) produces discon- 
tinuities in the phase characteristic associated with each 


crossing of the w axis by the resulting Fourier series.. 
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Fig. 4.5 Phase vs. normalized frequency. Non- 
recursive filter with coefficients obtained 
from Fourier series. 
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B. OYNTHESIS OF A SPECIAL TYPE OF RECURSIVE FILTER FROM 
A GiVEN MAGNITUDE SPECTRUM FUNGI=G) 


The method described in the previous section, where non- 
recursive filters were obtained for matching a given magnitude 
function, has been generalized by the author and allows a 
mecursive filter to be obtained. 

This special type of filter has a Z transform of the 


impulse response 


ata zt ae ee eee eee 2 entl ae ee 
H(z) n n-L 0 n-L < 
he be er ee een ém+l . 4 272m 
m-L oO m-L m 
CAs la 


AS CcanOmical form of this filter 1s “shown in’ Pigueewc. 

Note that the polynomials of the numerator and the 
denominator have an even number of delays for the input and 
output signals, and furthermore are mirror image polynomials 
because of the symmetry in the coefficients. Such a filter 


can be considered to be derived from the expression 





nueza B(z) (4.18) 
a + a2 Hyon Ey Fock By en. {cea 
where A(z) = 
(4.19) 
b + yz once Bae acct bau 2 ae +b, 2 
and 1) a 
(2220) 
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The frequency Specerum Of Such women sie coe erme,, 


7) Coe 
) = pceuams (4 ee) 


QJvT 
Eee) 


H ( 


emcee from Eqs. (4.12)and (4.13) 


=] a7 J (n-m) oT 
) = a ” (4.22) 


b + . 25, COS man r 

ie k 

Equation (4.22) shows the properties of this special type 
of filter being considered. It indicates that its frequency 
Spectrum is in the form of a ratio of two finite Fourier series, 
and the phase is linear with respect to frequency. 

When a magnitude versus frequency characteristic iS given 
to be approximated with a digital filter, if an approximation 
in the form of the ratio of two finite Fourier series can be 
obtained, then a method to synthesize the filter is available 
by using Eq. (4.22) as will now be explained. 


Hes The Expansion of a Maqnitude Function as a Ratio 
of Two Finite Fourier Series 


Consider first the properties of the Chebyshev 
polynomials, as discussed for instance by R. W. Hamming [13]. 
"The Chebyshev polynomials have all the properties 
of both the Fourier series and the orthogonal polynomials; they 
are, in fact, the Fourier functions cos n 6 in the disguise 


of a simple transformation of the variable 


6 = arc cos x 
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If a variable x 1s defined such that 


X = cos wT (4,23) 
iam waere =) < X < I 
then, 

TL) = cos (n arc cos x) = cos (wTn) (A424) 


In view of Eq. (4.24) an expansion in Chebyshev 
polynomials in the variable x corresponds to a Fourier series 
in the variable wT, when the two variables are related as in 
Eq. (4.23). This property can be used to approximate a given 
magnitude spectrum function as a ratio of two finite Fourier 
series in the following manner: 
oa 


Given a function |[G(e which has to be “appeerar. 


mated consider the following transformation 


\eiel™ 


)| = G, (wT) = G, (are cos x) = H(x) (4 2p 
If this is done numerically, that is if pairs of 
values describing points of the magnitude versus frequency 
Gharacteristic are given Eq. (4.25) can be considered as a 
change in the wT values of the given points, maintaining the 
ordinate values. The interest in doing this is that now the 
new function of x can be approximated as a ratio of two 
polynomials in x. Later the polynomials may be expressed as 
Chebyshev polynomials and then use of the property described 


om) (4.24) leads to the Fourier Series coefficients. —=raeess 
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Me Goee eam be approximated by 


ees, 
H (x) ~ O (x) (eZ 6) 


m 


and numerator and denominator polynomials can be represented 


as an expansion of Chebyshev polynomials of the form 


Ge) (4.27) 


Wsing the property of Eq. aaa) 


ec  — ——————— == (4.28) 


and Eq. (4.28) represents the approximation of a given function 
es the ratio of two finite Fourier series. 
For synthesis purposes, the coefficients in Eq. (4.28) 


can be equated with those of Eq. (4.22) yielding 
k # 0 (4.29) 
kes e0 (4.30) 


The complete synthesis procedure is illustrated in 


Fig. 4.9. The following points should be noted: 
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a. If the discrete steps in the abcissae of 
Ig (et) | are uniformly spaced, the abcissae of H(x) will be 
nonuniformly spaced since x) = cos OT. 

be Expressing H(x) = P(x)/Q(x) is possible in many 
ways. For the work presented here {BM subroutine ARAT is used. 

ore Expressing P(x) and Q(x) as Chebyshev poly- 
nomials is performed using the following tables taken from 


Mecracken and Dorn [14]. 


a Chebyshev Polynomials 


Ty (x) — dh 
T, (2) = 
- oo 
T., (x) Dae 1 
3 
T(x) — 4 oe eee 
T,(x) = 8 | = 760 a 
T. (x) = eG x - 20 re a sc 
Oe Powers of x as functions of TS () 
1 = To 
x = Ty 
7 ate 
x = 2 (T)+T.) 
Boas 
x = 4 (3T,+T.) 
a 
x = g (3T)+4T,+T,) 
5 1 
x = Te (10T, +5T,+T.) 


The degree of the numerator and denominator 


polynomials can be different, as indicated in Eq. (4.30) 2 Le 
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difference in polynomial degree has direct influence on the 
phase characteristics, as shown in@iq. (4522). (Selleetimi ee em 
degrees for numerator and denominator, 10 1s posstelte = rommae 
the phase equal to zero for all frequencies. Otherwise the 
phase will change linearly with frequency, and the slope will 
be directly proportional to the difference in degree.. 

An example of synthesSis uSing this procedure is now 
presented. For comparison with the M. Martin method, the 
example given in the previous section is reworked. The 
Magnitude spectrum of Fig. 4.3 will be approximated using a 
filter with eight delays, four in the numerator and four in the 
denominator, which implies the approximation of the given 
magnitude spectrum with a ratio of two Fourier series involving 
up to the second harmonic in each. The method used in the 
example to obtain the ratio of two Fourier series approximating 
fie given magnitude spectrum function, is the numerical method 
illustrated in Fig. 4.9, uSing subroutine ARAT of the IBM 
scientific package to obtain the ratio of two polynomials in 
the transformed variable x = cos wT. The results are: 

Given eight points aceeriee the magnitude spectrum 


function equally spaced on the wT axis in the range 
Tie SSS ay (A. 31) 


The change to the variable x is performed and a ratio of 
two polynomials in x approximating the points was found 


2 
ieee Ope Ola eee 105 UALS S4 ese ern ett x (4.32) 


1000 = 1.353494 x + "0247 720m 
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The polynomials are now expanded in terms of 
Chebyshev polynomials using table shown previously. This yields 
0301823 Ty (x) + O2,0241¢ T, (x) + 000 G20 T. (x) 


Bl) es (4.33) 
ee 3370 Ty (x) = des94g4 Ty (x) + O2250 70 T(x) 





Using Eq. (4.24) the ratio of Fourier series can be 


obtained. 
Ig (dT) | = 2:01823 + 0-02418 cos wT + 0.00639 cos 2uT 
1.238/0 -—- 1.39494 cos Gf + 0.233 /6 ces Jan 
(4.34) 
From Bqs. (4.22) and (4.28) 2tetollows that 
-1} -2 = ay, 
Giz) = 02003194+0.012092 7" +0.01823z “+0.012092 "+0.003192 
: : j 


0.11935-0.69747z “41. 239702 1-20 6o74 77 | Oeo en 


(4.35) 


The magnitude and phase of the spectrum has been 
mibotted and illustrated in Figs. 4.10a and 4.10b. As a@ comparrson 
with M. Martin's method shown in Fig. 4.4, it can be noted that 
there is an increase in attenuation and elimination of lobes 
as well as zero phase shift. 

Cc. SYNTHESIS OF RECURVSIVE FILTERS FROM MAGNTIDUE SQUAKED 
PReOURNGY SPECTRUM FUNCTIONS 
“In the previous section the filter is restricted to a 
minor image coefficient realization starting with the magnitude 
of a transfer function. In this section a new approach is 
developed to synthesize recursive digital filters, from a given 


Magnitude squared function. 
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Fig. 4.10b Phase vs. normalized frequency. Recursive 
filter, coefficients obtained from numerrteal 
approximation described in Fig. 4.8. 
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It 1s shown that a digital filter SGescribed anes -cene se 
general form has a magnitude squared frequency spectrum which 
is a ratio of two finite Fourier series with cosine terms only. 
A method for obtaining an approximation to such a function, 
has been described in the previous section, using a change of 
variable and expansion of the approximating function ina 
series of Chebyshev polynomials. In this section the coeffi- 
cients of the filter are obtained by equating the coefficients 
of the Fourier series expansion witn the coefficients of the 
filter according to some nonlinear relationships which are 
developed. 

ia Bove lopencierOr enue Lneom, 

Consider a generalized description of a digital 


fmmelbcer, OoL the form 





k : 
- _ oe mee 
ao ay a, 2 : + ay Z d +...+ a, 2 HE Be oar es 
H(z) = ———-—_ ———_._ -_—_————-- >: = ees (4.36) 
-l -2 -m™m -1 
b. + b. Zz + b. Z +...+ b Z > \opaeya 
0 1 2 m 0 ike 


The spectrum of this filter can be obtained by 
evaluating the magnitude and phase of Eq. (4.17) for values of 


measong the unit circle 


Z=e = cos wT + j sin wT CA ore 
Substituting (4.37) inte ssG). yierds 


eeclo S as aoe IP (ek = oe 


| * - (4.38) 
go oe Eaae ee eer 
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Consider 


Jot _ N(Gje) 
Her) = Da) 


; Z 


|D(jw) | 7 


. 2 Z 
a ee 
Hew ee eo AE 7 AL (4.39) 
Re {D(jw) } + Im” {D(jw) } 
ns 2. s 2 
(at cos nwT) + (2 ae Sin nwT) 
Ju(eJ7) |* = ———et net (4.40) 
(b +2 cos nee + (2 b sin nw?) * 
‘ai — IP n=1 ” 
which can be reduced to the form 
k 5 k~-1 k-2 
> a +2} aa +1 cos wWT+2t aa +2 cos 2wT +... 
jwT, 2 _ n=0 n=0 n=0 ” 7 
|H(e i = m m-l —2 7 
5 m 
® jb “+2% bb 1 cos wT+2r bb cos 2wT +... 
watt n 5 n nt iy n n+2 
leet) 
k 
| " GC Cos nul 
Ju(e@I°F) |* = BOP __ (4.42) 
2; ai ees not 
n 
n=0 


Proof of Eq. (4.41) is given in Appendix A.. Equation 
(4.41) indicates that the magnitude squared frequency spectrum 
is the ratio of two finite Fourier series with cosine terms 
only. Equation (4.42) shows the relationships between the 
coefficients of the Fourier series, om and do and the coeffi- 


cients of the filter ee and b.. Thus, from (4.22) and (4723) 
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Mm 
Cn = 2 : oe Den (4.435b) 


This relationship can be written in matrix form as 


follows 
Ay 4, Ay seeeee ay. ay Co 
0 ay Ay seeeee Ay ay ; C12 
0 0 Ay reece. - AL _o a. £272 (4.44) 
0 0 CS Mice ge Ay a, | bey 3 
Similar relationships hold for the denominator 
coefficients. As discussed in the previous section, it is 


possible to find an approximation to an even function as a 
rational expression of Fourier series so that Ch and qd are 
known. The coefficients of the Fourier series, and hence the 
vector of c's in Eq. (4.44), will be obtained from the expansion 
of the desired spectrum function. Solving for the a's in 
me. (4.44) yields the coefficients of the numerator of the 
Gigital filter realization. An identical procedure can be 
followed for the coefficients of the denominator. 

This set of nonlinear equations has several solutions, 
in fact, as k gets larger, the number of solutions also 


increases. This is not surprising because of the nature of the 
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problem. That is, given a filter, a unique magnitude squared 
frequency spectrum function can be found, but the inverse 
procedure is not unique. Given the magnitude of the spectrum, 
the descriotion of the filter is incomplete because there is 
no information about the phase characteristics. 

Consider an example of the use of this synthesis 
method to approximate a given magnitude squared function. For 
comparison with the example given the same magnitude spectrum 
Shown in Fig. 4.3 will be squared. What results in the same 
function, which is now approximated with a filter using only 
two delays for the output and input signals. The nonlinear 


equations to be solved are of the form 





2 2 a ; 
ay + ay = as = Co G4 .45) 
ay ay + ay a. = C12 (4.46) 
a, a = Co/2 (die a7 } 
a, = (4.48) 
; Cc 
ne ¢ - °2 (4..49) 
4a, 4a, 2 
Cc 
— _- (4.50) 
2 2a 
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Equation (4.48) indicates there are four solutions 
£Oxr ay and Eq. (4.49) indicates there are two solutions of Ag 


for each ay: Hence there are eight sets of coefficients 
Satisfying the given nonlinear equations.) [he \same etc meare 
for the denominator set of coefficients and therefore there 
are 64 solutions for the filter. In general not all the 
solution may be valid if only real coefficients are being 
considered. 

Using the results obtained in previous example, 
Eq. (4.34), the magnitude squared function is expanded in the 
same ratio of Fourier series, yielding 
jwT iC 0.01623 +.0.02418 cops eee OCG 3S rcosc Zak 


ne ee «6(4,51) 


| H(e 1.23870 — 1.39494 cos @& $90 253/0 cos Zur 


) 


Substituting these coefficients into (4.48), (4.49), 
(250), two sets of solutions are obtained, called filter A 


and B respectively. 


FILTER A FILTER B 
ay = 0.038983093 0.081958605 
a, = 0.099965521 0.099965521 
a, = 0.081958605 0.038983093 
b, = -0.284898296 -0.418921424 
b, = 0.990978207 0.990978207 
b, = -0.418921424 -0.284898296 


The magnitude and phase versus frequency characteristics 


ome shownein Pag. 4:11, 4.12, 4.13, 4.14. 
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magnitude squared criterion. Filter B. 
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Comparison of filters A and B indicates the same 
Magnitude for both filters and a different phase characteristic 
aS was expected. With these sets of coefficients, still two 
more combinations are possible, that is the numerator of A 
and the denominator of B and vice versa. At this point, 
nothing has been said about the phase characteristic of this 
synthesis method. It is left as a suggestion for further 


research, the introduction of phase characteristics restrictions. 
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V. SUMMARY AND SUGGESTIONS 
FOR FURTHER RESEARCH 

The two main new ideas introduced in the thesis are: 

i A generalization of Martin's synthesis technique to 
include recursive filters, which allows filters with zero 
phase or linear phase characteristics, to be oktained. 

2 A new synthesis technique, which depends upon the 
expansion of a given function as a ratio of two Fourier series. 
The resulting technique involves the solution of simultaneous 
nonlinear algebraic equations, but the filter requires half 
the storage of Martin's method for the same degree of approxi- 
mation. The number of multiplications to be performed is also 
decreased by one half. 

The foregoing studies have given rise to several basic 
euestions for further research: 

ANA Explore the concept of expanding a given function 
as a ratio of two finite Fourier series directly from Fourier 
series expansion theory, as versus the change of variable 
technique, involving Chebyshev polynomials, developed in this 
thesis. It would appear that an unlimited number of ratios 
are possible. If this is the case, filter stability can be 
assured by selecting the poles of the digital transfer function 
to lie inside the unit circle of the z plane, and then selecting 
the coefficients of the numerator to obtain the desired 


magnitude spectrum. 
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Dee exionder synthesis procedures starting from a class- 
ical specification for the filter in terms of tolerance bands 
and attenuation rates. Determine how these design criteria 
ean be met. 

oe Find a general method to solve the nonlinear equa- 
tions, Eq. (4.44), which were developed for the filter synthesis 
proceedure starting from a magnitude squared function as 


aascussed. 
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Appendix A 


Derivation of Equation (4.22) 


‘oT. 12 k =» Js , 
lu (e? ) | = (ay + 2% a cos nuwT)~+(2z Sin nwT) (A.1) 
n 
= n=l 
k > k 5 
bo + 2 bo cos nuwT) ~+(2 Sin mwT) 
=] n=] ™ 
consider the numerator only 
2 us on 8 | 2 
IN(wT)]“ = a, + 2 a. cos nwT)“ + (Z a. sin nwT) 
0 n n 
2 i s 2 us 2 
= a. +2a,% a cos nwT+(Z a cos nuwT) +(2 awe sin nwT) 
0 Ce - se 
n=l n=l n=l 
(A.2) 
The last two terms can be expanded as follows: 
: A | Ro eg 
(x a, cos nwT)" + (2 a, Sin nwT)” = 2 a, cos nwT 
=] n=1 n=] 
ea k k 5 
+ 2y a_cos nuT 2 a. COSmaiaiedt. 2 a Ssin® mut 
=] j=nt+] n=] 
k-1 k 
oZs ea sin nul a. Sin jwT (A. 3) 
=1 j=nt+l 
2 eo k 
= Bera. (cos nwoltsin wy mol) 2 eeecces nu) a. cos jwT 
n 
n=l n=l j=entlL 
s . 
- cin nul a. sin jwT } (A. 4) 
n=n+l 
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k 5 dl as 
= 3 a, + oe {2 a, cos nwT a. cos jwT 
n=1 n=l * j=n+1l J 
k 
+ a cos nwT a. sin juT (A.5) 
n=n+1 J < 
k 5 kes 
= ya ae » aa. {cos nwT cos juT 
n=1 n=l ‘ j=n+l 
+ sin nwT sin jot} + (A.6) 


Use the trigonometric identity 


{cos nwT cos jwT + sin nwT sin jwT} = cos (j-n)wT Oe) 
k 5 Ik IR a 
= 3 a, t+ oe {2 aa, Cos Se 
i=l n=l ‘j=nt+1 ” J 
eianging subscripts, 1 = j-n , j = n¥i the limits of the 


second summation become 


ee peal ig i=l 
J =k 2 1 = k-n 
k 5 k-l k-n 
= > a + 2 f n a a ,; Cos Teh 6: Warsp 
n=l 2 hee or 


Changing summation order Eq. (A.8} becomes 


k 5 eth Sah : 
= Lea +2 F kK ana _ } cos Toto (A.9) 
n ; n nti 
n= | Te 


or 





Combining (A.9) with (A.2) yields 


2 2 i | k 2 
[N(wr) |“ = ay +22  apap,, cosiwt +2 a 
i=] n=] ” 
k-l1 k-l 
+ 2 a be ananeg Cos 1 wT (Aso 
i=l n=l 
k 3 k-1 
= -_ + 2 Ay ay cos k wT + 2 - Aya ouG COS Gar 
k-l1 k-l 
+ 2 °F , a.a@ .. cos i u0F (A.11) 
_, n nti 
i-l n=l 
k 5 k-1l1 k-i 
= ya + 2 a,a, cos k wT + 2 Y a.a.. Coo 
= © sigs i=l n=o 7 ** 
(Awl 2) 
Ne 5 Ke k-1 
= =) oe % a a... 0s 2.7002 (A.13) 
oe ia _ 7 n nti 
n=0 i=l n=0 


The same relationship holds for the denominator, and so 


established the validity of Eq. (4.22). 
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